scRNAseq long read analysis - Part 2

In Part 1, we showed how to preprocess scRNA-seq isoform datasets and create a multi Assay Seurat object. In this part, we’ll use this Seurat object to perform differential isoform expression analysis and explore pertinent biological questions. We’ll focus on searching for genes with a transcript isoform expression that differed between clusters, which we refer to as “isoform switches”.

We start the analysis by loading the necessary packages:

library(dplyr)
library(patchwork)
library(ggplot2)
library(reactable)

Plotting functions

see_isoforms_expression = function(sobj, gene, limits) {
  isoforms = grep(paste0("^", gene, "\\."), rownames(sobj), value = TRUE)
  
  plot_list = lapply(isoforms, FUN = function(one_isoform) {
    p = Seurat::FeaturePlot(sobj, features = one_isoform,
                            reduction = name2D, order = TRUE) +
      ggplot2::labs(title = one_isoform) +
      ggplot2::scale_color_gradientn(limits = limits, colors = c("lightgray", "#FDBB84", "#EF6548", "#7F0000", "black")) +
      ggplot2::theme(aspect.ratio = 1,
                     plot.subtitle = element_text(hjust = 0.5)) +
      Seurat::NoAxes()
    return(p)
  })
  
  return(plot_list)
}

see_genes_expression <- function(gene){
  Seurat::FeaturePlot(seurat_obj, gene, reduction = name2D)+
      ggplot2::scale_color_gradientn(limits = c(0,4), colors = c("lightgray", "#FDBB84", "#EF6548", "#7F0000", "black")) +
      ggplot2::theme(aspect.ratio = 1,
                     plot.subtitle = element_text(hjust = 0.5)) +
      Seurat::NoAxes()
}

We set the global settings of the analysis. We will store data there :

out_dir = "./output/"
data_dir <- "./data/"

2.1 Loading the Seurat object

We load the Seurat object:

seurat_obj = readRDS(paste0(out_dir, "/datasets/", "seurat_obj_annotated.rds"))
color_markers = readRDS(paste0(out_dir, "/datasets/", "color_markers.rds"))

This is the projection of interest:

name2D =  "RNA_pca_10_tsne"

2.2 Cell visualization

We visualize cells:

cell_type = Seurat::DimPlot(seurat_obj, reduction = name2D, group.by = "cluster_cell_type") +
  Seurat::NoAxes() +
  ggplot2::scale_color_manual(values = unlist(color_markers),
                              breaks = names(color_markers)) +
  ggplot2::theme(aspect.ratio = 1,
                 plot.title = element_text(hjust = 0.5),
                 plot.subtitle = element_text(hjust = 0.5))

clusters = Seurat::DimPlot(seurat_obj, group.by = "seurat_clusters", label = TRUE,
                           reduction = name2D) +
  Seurat::NoAxes() + ggplot2::ggtitle("Clustering") +
  ggplot2::theme(aspect.ratio = 1,
                 plot.title = element_text(hjust = 0.5))

clusters | cell_type

2.3 Differential isoforms expression (DIE)

In this section, we’ll use the ISO_SWITCH_ALL function from Isoswitch to analyze isoform switches between cell types. we refer to “isoform switches.” when two isoforms of the same gene are markers of different cell clusters. The ISO_SWITCH_ALL function uses Seurat’s FindMarkers function to identify marker isoforms for each cluster in a gene. It outputs a data frame, listing each transcript marker in a separate row.

clusters = levels(seurat_obj@active.ident)
switch_markers = isoswitch::ISO_SWITCH_ALL(seurat_obj, clusters, assay = "ISO", 
                                           min.pct = 0, logfc.threshold = 0.40)

2.3.1 Save DIE list

We save the list of DE isoforms :

saveRDS(switch_markers, file = paste0(out_dir, "./datasets/", "switch_markers_combined_all.rds"))

2.4 DIE visualisation

We load the DE isoforms table :

switch_markers = readRDS(paste0(out_dir, "/datasets/switch_markers_combined_all.rds"))

2.4.1 Explore the switch markers table

In this section, we explore the switch_marker dataframe.

List of isoform switchs:

isoswitch::gene_switch_table(switch_markers)

Can you identify the genes and cell types that show the most significant changes in isoform expression?

To facilitate the graphical interpretation of the results:

• plot_marker_matrix() produces a heatmap of number of unique genes per contrast between clusters
• plot_marker_score() produces a volcano plot showing p-values and average logFC for each gene with an isoform switch

We visualize the number of DE isoforms between each cluster :

pl1 = isoswitch::plot_marker_matrix(seurat_obj, switch_markers)
pl2 = isoswitch::plot_marker_score(seurat_obj, switch_markers)
pl1 | pl2

Individual volcano plots for each cluster.

isoswitch::plot_marker_score(seurat_obj, switch_markers, facet=TRUE, ncol=3)

2.4.2 Visualization of specific genes

In this section, we visualize isoforms expression level for specific genes.

First, we set the assay to “ISO” because all expression levels are stored there.

Seurat::DefaultAssay(seurat_obj) = "ISO"

Let’s take the example of Clathrin light chain A (Clta) gene isoform expression switch during neuronal maturation:

see_isoforms_expression(seurat_obj, "Clta", c(0,4)) %>%
  patchwork::wrap_plots(., ncol = 3)

Seurat::FeaturePlot(seurat_obj, reduction = name2D, features = c("Clta..ENSMUST00000107849", "Clta..ENSMUST00000170241"), blend = TRUE, order = TRUE)

2.4.3 Report

In this section, we generate a small report.

We load annotation files: The Isoswitch documentation provides a comprehensive guide on how to build your annotation metadata.

gtf_df <- readRDS(paste0(data_dir,"annotation/gtf_rds.rds"))
gene_metadata <- readRDS(paste0(data_dir,"annotation/gene_metadata.rds"))
transcript_metadata <- readRDS(paste0(data_dir,"annotation/transcript_metadata.rds"))

For Clta gene:

gene_switches <- isoswitch::compute_switches(switch_markers, gene="Clta")
knitr::kable(isoswitch::format_switch_table(gene_switches))

geneId	t1	c1	p1	log2fc1	t2	c2	p2	log2fc2
Clta	ENSMUST00000107849	Glutamatergic	2.69e-23	5.35	ENSMUST00000170241	radial_glia	4.11e-20	-2.29
Clta	ENSMUST00000107849	Glutamatergic	4.75e-27	2.35	ENSMUST00000170241	cyclin_radial_glia	7.11e-20	-2.03
Clta	ENSMUST00000107849	imature_Glutamatergic	1.39e-20	5.03	ENSMUST00000170241	radial_glia	3.87e-12	-1.55
Clta	ENSMUST00000107849	GABAergic	1.09e-18	5.54	ENSMUST00000170241	radial_glia	1.56e-11	-1.73
Clta	ENSMUST00000107849	GABAergic	2.35e-19	2.54	ENSMUST00000170241	cyclin_radial_glia	1.34e-10	-1.46
Clta	ENSMUST00000107849	imature_Glutamatergic	9.33e-20	2.03	ENSMUST00000170241	cyclin_radial_glia	1.54e-10	-1.28
Clta	ENSMUST00000170241	intermediate_progenitors	3.60e-10	-2.11	ENSMUST00000107849	Glutamatergic	4.35e-09	1.21
Clta	ENSMUST00000107849	GABAergic	3.37e-06	1.40	ENSMUST00000170241	intermediate_progenitors	1.55e-05	-1.54
Clta	ENSMUST00000107849	GABAergic	1.09e-18	5.54	ENSMUST00000107845	radial_glia	8.57e-05	-2.16
Clta	ENSMUST00000170241	intermediate_progenitors	1.50e-04	-1.36	ENSMUST00000107849	imature_Glutamatergic	8.12e-03	0.89
Clta	ENSMUST00000107849	imature_Glutamatergic	1.39e-20	5.03	ENSMUST00000107845	radial_glia	1.57e-04	-1.84
Clta	ENSMUST00000107849	Glutamatergic	2.69e-23	5.35	ENSMUST00000107845	radial_glia	4.01e-03	-1.94
Clta	ENSMUST00000107849	imature_GABAergic	5.20e-07	4.49	ENSMUST00000170241	radial_glia	7.68e-03	-0.81
Clta	ENSMUST00000107851	Cajal_Retzius	2.21e-02	6.09	ENSMUST00000170241	radial_glia	2.58e-02	-2.45
Clta	ENSMUST00000107849	Cajal_Retzius	2.43e-04	4.38	ENSMUST00000170241	radial_glia	2.58e-02	-2.45
Clta	ENSMUST00000107849	Glutamatergic	2.46e-09	-0.85	ENSMUST00000170241	imature_GABAergic	2.70e-02	1.49
Clta	ENSMUST00000107849	imature_GABAergic	5.20e-07	4.49	ENSMUST00000107845	radial_glia	2.97e-02	-1.38

isoswitch::isoswitch_report(seurat_obj, "ISO", gene="Clta", marker_list=switch_markers, gtf_df, transcript_metadata) 

## [1] "Locus"
## [1] "Junction"
## [1] "umi counts"
## [1] "dotplot"
## [1] "PATCHWORK"

• Now, try applying the same analysis to explore isoform expression for other genes of interest.

2.5 Save

We save the list of DE isoforms :

saveRDS(switch_markers, file = paste0(out_dir, "./datasets/", "switch_markers_combined_all.rds"))

2.6 Ressources

• Seurat tutorials
• Isoswitch R package
• Read more about Single-cell long-read

R session

## R version 4.2.2 (2022-10-31)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS 14.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] reactable_0.4.4 ggplot2_3.4.0   patchwork_1.1.2 dplyr_1.1.4    
## 
## loaded via a namespace (and not attached):
##   [1] backports_1.4.1        spam_2.10-0            Hmisc_4.7-2           
##   [4] plyr_1.8.8             igraph_2.0.1.1         lazyeval_0.2.2        
##   [7] sp_1.5-1               splines_4.2.2          RcppHNSW_0.5.0        
##  [10] crosstalk_1.2.0        listenv_0.9.0          scattermore_1.2       
##  [13] digest_0.6.31          htmltools_0.5.4        fansi_1.0.3           
##  [16] checkmate_2.1.0        magrittr_2.0.3         tensor_1.5            
##  [19] cluster_2.1.4          ROCR_1.0-11            limma_3.54.2          
##  [22] globals_0.16.2         matrixStats_0.63.0     spatstat.sparse_3.0-0 
##  [25] jpeg_0.1-10            colorspace_2.0-3       ggrepel_0.9.2         
##  [28] xfun_0.36              jsonlite_1.8.4         progressr_0.12.0      
##  [31] spatstat.data_3.0-0    survival_3.4-0         zoo_1.8-11            
##  [34] glue_1.6.2             polyclip_1.10-4        gtable_0.3.1          
##  [37] leiden_0.4.3           ggtranscript_0.99.9    future.apply_1.10.0   
##  [40] isoswitch_0.0.0.9000   abind_1.4-5            scales_1.2.1          
##  [43] DBI_1.1.3              spatstat.random_3.0-1  miniUI_0.1.1.1        
##  [46] Rcpp_1.0.9             viridisLite_0.4.2      xtable_1.8-4          
##  [49] htmlTable_2.4.1        reticulate_1.34.0      foreign_0.8-84        
##  [52] dotCall64_1.1-0        Formula_1.2-4          htmlwidgets_1.6.0     
##  [55] httr_1.4.4             RColorBrewer_1.1-3     ellipsis_0.3.2        
##  [58] Seurat_5.0.0           ica_1.0-3              pkgconfig_2.0.3       
##  [61] farver_2.1.1           nnet_7.3-18            sass_0.4.4            
##  [64] uwot_0.1.14            deldir_1.0-6           utf8_1.2.2            
##  [67] tidyselect_1.2.0       labeling_0.4.2         rlang_1.1.3           
##  [70] reshape2_1.4.4         later_1.3.0            reactR_0.5.0          
##  [73] munsell_0.5.0          tools_4.2.2            cachem_1.0.6          
##  [76] cli_3.6.2              generics_0.1.3         ggridges_0.5.4        
##  [79] evaluate_0.19          stringr_1.5.0          fastmap_1.1.0         
##  [82] yaml_2.3.6             goftest_1.2-3          knitr_1.41            
##  [85] fitdistrplus_1.1-8     purrr_1.0.2            RANN_2.6.1            
##  [88] pbapply_1.6-0          future_1.30.0          nlme_3.1-161          
##  [91] mime_0.12              compiler_4.2.2         rstudioapi_0.14       
##  [94] plotly_4.10.1          png_0.1-8              spatstat.utils_3.0-4  
##  [97] tibble_3.2.1           bslib_0.4.2            stringi_1.7.8         
## [100] highr_0.10             RSpectra_0.16-1        forcats_0.5.2         
## [103] lattice_0.20-45        Matrix_1.6-2           vctrs_0.6.5           
## [106] pillar_1.9.0           lifecycle_1.0.3        spatstat.geom_3.0-3   
## [109] lmtest_0.9-40          jquerylib_0.1.4        RcppAnnoy_0.0.20      
## [112] data.table_1.14.6      cowplot_1.1.1          irlba_2.3.5.1         
## [115] httpuv_1.6.7           R6_2.5.1               latticeExtra_0.6-30   
## [118] promises_1.2.0.1       KernSmooth_2.23-20     gridExtra_2.3         
## [121] parallelly_1.33.0      codetools_0.2-18       assertthat_0.2.1      
## [124] fastDummies_1.7.3      MASS_7.3-58.1          withr_2.5.0           
## [127] SeuratObject_5.0.0     sctransform_0.4.1      parallel_4.2.2        
## [130] reactablefmtr_2.0.0    grid_4.2.2             rpart_4.1.19          
## [133] tidyr_1.3.1            rmarkdown_2.19         Rtsne_0.16            
## [136] spatstat.explore_3.0-5 shiny_1.7.4            base64enc_0.1-3       
## [139] interp_1.1-3